Cite this Article
Wang Jing-Jing, Yao Yang, Gao Zhi-Wei, Li Xiao-Li, Wang Jun-Song. Generating mechanism of pathological beta oscillations in STN–GPe circuit model: A bifurcation study. Chinese Physics B, 2020, 29(5): 058701  
Permissions
Generating mechanism of pathological beta oscillations in STN–GPe circuit model: A bifurcation study
Wang Jing-Jing1, Yao Yang1, Gao Zhi-Wei2, Li Xiao-Li3, Wang Jun-Song1, †
School of Biomedical Engineering, Tianjin Medical University, Tianjin 300070, China
Faculty of Engineering and Environment, Northumbria University, NE1 8ST, UK
National Key Laboratory of Cognitive Neuroscience and Learning, Beijing Normal University, Beijing 100875, China

 

† Corresponding author. E-mail: wjsong2004@126.com

Project supported by the National Natural Science Foundation of China (Grant Nos. 61473208 and 61876132) and the Tianjin Research Program of Application Foundation and Advanced Technology, China (Grant No. 15JCYBJC47700).

Abstract

Parkinson’s disease (PD) is characterized by pathological spontaneous beta oscillations (13 Hz–35 Hz) often observed in basal ganglia (BG) composed of subthalamic nucleus (STN) and globus pallidus (GPe) populations. From the viewpoint of dynamics, the spontaneous oscillations are related to limit cycle oscillations in a nonlinear system; here we employ the bifurcation analysis method to elucidate the generating mechanism of the pathological spontaneous beta oscillations underlined by coupling strengths and intrinsic properties of the STN–GPe circuit model. The results reveal that the increase of inter-coupling strength between STN and GPe populations induces the beta oscillations to be generated spontaneously, and causes the oscillation frequency to decrease. However, the increase of intra-coupling (self-feedback) strength of GPe can prevent the model from generating the oscillations, and dramatically increase the oscillation frequency. We further provide a theoretical explanation for the role played by the inter-coupling strength of GPe population in the generation and regulation of the oscillations. Furthermore, our study reveals that the intra-coupling strength of the GPe population provides a switching mechanism on the generation of the abnormal beta oscillations: for small value of the intra-coupling strength, STN population plays a dominant role in inducing the beta oscillations; while for its large value, the GPe population mainly determines the generation of this oscillation.

PACS: ;87.19.lj;;87.19.ln;;82.40.Bj;
Keyword:;STN-GPe model;;coupling strength;;beta oscillation;;bifurcation;
1. Introduction

Parkinson’s disease (PD) is a common chronic neurodegenerative disorder that mainly occurs in the elderly and has an adverse effect on the quality of the patient’s life. The PD is attributed to the degeneration and death of dopaminergic neurons in the midbrain and the substantia nigra pars compacta.[1] Such a pathological degeneration results in the decrease of the dopamine (DA) content of the brain, which causes the abnormal discharges to occur in BG nuclei of both patients[2] and animal models[3] of PD. Hence, the decrease at the level of dopamine severely affects the neuronal activities in the cortico-basal ganglia (BG)-thalamocortical loops.[4,5] As is well known, the rigidity and bradykinesia[6] are the main Parkinsonian symptoms, which are closely correlated with oscillations in the beta band frequency (12 Hz–35 Hz) in the basal ganglia.[7]

The exact pathophysiological origin for the appearance of the synchronized beta oscillations remains controversial.[8–14] Some studies have shown that the pathological beta oscillations were generated in the striatum and then spread to the rest of the loop.[12] The striatal theory suggested that the STN–GPe model could generate beta frequency oscillations when the input from striatum to GPe increased.[9,11] Others suggested that beta oscillations originated from the cortico-thalamic circuit and spread to the BG nuclei,[8] and there were also studies suggesting the beta oscillations were generated internally within the BG network model.[10,13,14] Overall, most of previous studies have shown that beta oscillations can be observed in the dopamine-depleted BG nuclei composed of STN and GPe, and in the striatum and cortex.[15,16]

Computational research is becoming a promising method to elucidate the mechanism behind the various neurological diseases, such as Parkinson’s disease[5,17,18] and absence seizures.[19–21] The STN and GPe in the basal ganglia connect to each other to construct a negative feedback loop, a typical neural motif for the generation of limit cycle oscillations. Thus, the STN–GPe feedback circuit model describes the possible origin of the abnormal beta rhythms.[10,13,14,22] The study of Holgado et al.[23] indicated that the coupling strengths between STN and GPe need to be sufficiently strong for the generation of beta oscillations in the STN–GPe model. Fei Liu et al.[24] proposed a neuron mass model of the basal ganglia composed of STN and GPe networks to reproduce spontaneous beta oscillations related to Parkinson’s disease by changing the coupling strengths and intrinsic parameters of the model.

Spontaneous oscillations are a special behavior in nonlinear systems, and related to limit cycle oscillations, an intrinsic phenomena of nonlinear dynamical systems.[25,26] The generation of spontaneous oscillations in the STN–GPe circuit model is primarily attributed to the sigmoid function.[27] In the study of Holgado et al.,[23] for the simplicity of analysis, the sigmoid functions in both STN and GPe populations are linearized with a linear function. Thus, the linear stability theory is used to derive the conditions for the generation of beta oscillations in the subthalamic nucleus–globus pallidus network. In the study of Fei Liu et al.,[24] to use the describing function method to explain the generation of excessive beta oscillations, the sigmoid function in the GPe population is approximated by its equivalent gain. However, replacing the sigmoid function by its linear form[23] or equivalent gain[24] may result in analysis errors due to the model approximating error. Moreover, most of previous studies conducted simulations to determine the boundary between the regions generating normal and abnormal beta activities in the two-dimensional model parameter space.[23,24] Thus, it still remains to be addressed by using a nonlinear analysis method to probe the dynamical mechanism behind the spontaneous beta oscillations in the STN–GPe circuit model.

In addition, Liu et al. found that the upper beta frequency originates from the high frequency cortical input, and there is a transition mechanism between the upper and lower beta oscillatory activities underlined by inhibitory self-feedback within the GPe.[27] Furthermore, the study of Liu et al. indicated that the STN primarily influenced the generation of the beta oscillation while the GPe mainly determined its frequency.[24] It is still unknown why the two populations play different roles in generating the abnormal beta oscillations, and what is responsible for the transition mechanism between the upper and lower beta oscillations.

A bifurcation analysis method, as an effective means to research the dynamic behavior of nonlinear systems, has been extensively used to study the neurological diseases,[25,28–30] such as epileptic seizures.[29,31] Here we employ the bifurcation analysis method to elucidate the mechanism for generating the pathological spontaneous beta oscillations underlined by coupling strengths and intrinsic properties of the STN–GPe circuit model. Moreover, by plotting frequency distribution in the two-dimensional model parameter space, we further reveal the regulating mechanism of the spontaneous oscillations frequency.

The rest of this paper is organized as follows. In Section 2, the STN–GPe feedback circuit model is described, which can reproduce pathological beta oscillations within certain parameter regions. In Section 3, by combining bifurcation analysis results and frequency distribution diagrams, we explore the influences of coupling strengths, synaptic and intrinsic parameters on the generation and regulation of the spontaneous beta oscillations in the model. We especially focus on the role of intra-coupling (inhibitory self-feedback) strength of GPe population played in generating the abnormal beta oscillations. Finally, the discussion and conclusion are given in Section 4.

2. Model and methods
2.1. STN–GPe circuit model

The STN–GPe feedback neural circuit model is proposed in previous studies to investigate the generation of pathological beta oscillations related to Parkinson’s disease (PD).[23,24,27] In the present study, this model is employed to further investigate the model parameters underlined generation mechanisms of beta oscillations associated with Parkinson’s disease. The architecture of the STN–GPe model is given schematically in Fig. 1(a). As shown in Fig. 1, the STN–GPe circuit model is mainly composed of excitatory STN and inhibitory GPe populations. In the model, STN projects excitatory connections to GPe, and receives inhibitory feedback connections from GPe with inhibitory self-feedback connections as well. The parameter C1 describes the inhibitory coupling strength from GPe to STN, and the excitatory coupling strength from STN to GPe is represented by parameter C2. The parameter C3 represents the self-feedback (intra-coupling) strength within the GPe. ySTN(t) is the output of the model, which is explained as a local field potential (LFP) signal. The Pg is the inhibitory external input of the GPe, and Ps represents the excitatory external input of STN. Here, we mainly study the generating mechanism of abnormal beta oscillations underlined by the model itself. Thus we set the values of both Ps and Pg to be zero.

Fig. 1. (a) Schematic diagram of STN–GPe circuit model. Red lines with rhombus represent the excitatory projections and black lines with round heads denote inhibitory projections. (b) Block diagram of STN–GPe circuit model, where C1 and C2 are the inter-coupling strengths between STN and GPe populations, and C3 represents the intra-coupling (self-feedback) strength of GPe population, and ySTN and yGPe are the outputs of STN and GPe populations, respectively. Both STN and GPe populations are denoted by their linear convolution operators (TSTN and TGPe) and sigmoid nonlinear elements (SSTN and SGPe).

Both STN and GPe neuronal populations are composed of the linear convolution operator TSTN(t) or TGPe(t) modeling excitatory or inhibitory synapse, and the sigmoid nonlinear function modeling spike generation close to the populations. The mathematical formulation at the population level is extensively used in some previously proposed models,[24–27,32–34] which were originally used to describe the neural activities of cortex.

The mean firing rate of every population is computed from its average membrane potential by the sigmoid function. Thus, for STN and GPe, the two functions of sigmoid shape are given by

where SSTN(v) and SGPe(v) are the firing states of the STN and GPe populations; v represents the input of the sigmoid element; rs and rg are the slopes of the sigmoid function at the origin representing the impacts of dopamine on STN and GPe, respectively,[35] es and eg are the maximum firing rates of the sigmoid functions of the STN and GPe populations, respectively.

For both STN and GPe populations, the synapses transfer mean firing rates of every population into their average membrane potential modeled by the two linear convolution operators TSTN(t) and TGPe(t), in the form of second-order linear functions as follows:

where Hs and Hg represent the maximum amplitudes of STN and GPe synaptic gains, and τs and τg are the time constants of passive membranes in STN and GPe, respectively. Thus the transfer functions of STN and GPe populations in the frequency domain are in the following forms:

Given this architecture, there are two main variables in the model: the outputs of the two synaptic boxes TSTN and TGPe in Fig. 1, denoted as y1(t) and y2(t), respectively, corresponding to the average membrane potentials of the two synapses in the STN and GPe populations. These variables satisfy the following differential equations:

All the parameters of the STN–GPe circuit model are given in Table 1.[24,27,36–40]

Table 1.

Parameter interpretation and values in STN–GPe feedback circuit model.

.
2.2. Methods

Bifurcation analysis is an efficient nonlinear method to elucidate the regulating mechanisms of the dynamics of a nonlinear dynamical system, and extensively employed to study the effects of model parameters on the dynamics of various neural models.[28–30] In the present study, to detect the beta pathological oscillations in the STN–GPe model, we conducted bifurcation analysis with respect to parameters of the model by using XPPAUT, a valid tool for simulating and analyzing dynamical systems.[41] To explore the influence of model parameters on the oscillation frequency, frequency curves with respect to model parameters and frequency distribution diagrams in two-dimensional parameter plane were plotted by using the MATLAB software.[42] The Runge–Kutta numerical integration method could be used to solve the differential equations at a zero initial condition.

To probe the role played by STN and GPe populations in generating the limit cycle oscillation, we define a measure named effect factor K to quantitatively evaluate the effect of the model parameters on the generation of spontaneous oscillations. The effect factor is defined as K = Δ2/Δ1, where Δ1 represents the change of the model parameter value, and Δ2 the change of the spontaneous oscillation region caused by Δ1. In Fig. 2, we suppose that the black line represents the stability boundary in two-dimensional plane (x,y), and the region to the right of the boundary line is the limit cycle regions, where x and y are the model parameters. Supposing the model parameter x changing from x1 to x2, i.e., Δ x = x2 − x1, which results in the fact that the oscillation region changes from y1 to y2, for example, Δy = y2 − y1, such that the effect factor of the parameter x on the generation of the limit cycle oscillation is formulated as Kx = Δy/Δ x = L, where L is the slope of the stability boundary. Obviously, the effect factor of the parameter y on the generation of spontaneous oscillation is formulated as Ky = Δx/Δ y = 1/L.

Fig. 2. Definition of effect factor K of model parameters on generating spontaneous oscillation.
3. Results

We first conduct codimension-one bifurcation analysis to study the influence of coupling strengths between STN and GPe populations on the spontaneous beta oscillations. In addition, through codimension-two bifurcation analysis, we explore how the model parameters interact to exert an influence on the generation of abnormal spontaneous beta oscillations, including inter-coupling strengths (for example, C1 and C2) between STN and GPe populations, and the properties of STN and GPe populations itself, such as synaptic gains (i.e., Hs and Hg) and sigmoid elements (for example, es and eg, and rs and rg). Especially, we focus on exploring how the self-feedback of the GPe (i.e., C3) exercises an influence on the spontaneous beta oscillations.

3.1. Influence of coupling strengths on pathological beta oscillations

Previous experimental and simulation studies have suggested that inter-coupling connection between STN and GPe and inhibitory intra-coupling connection in GPe play important roles in inducing beta frequency band behavior in the STN–GPe model.[23,24] In this subsection, to explore the influence of the coupling strengths on the generation of pathological beta oscillations, we conducted bifurcation analysis and plotted the curves of the frequency versus the coupling strength, while the other parameters had default values as listed in Table 1.

3.1.1. Influence of coupling strength on pathological beta oscillations

We first plotted the codimension-one bifurcation diagram versus coupling strength C, letting C2 = C, C1 = 0.5C, and C3 = 0.5C, which means that the inter-coupling strength between STN and GPe populations and intra-coupling strength of GPe population change in the same way. The bifurcation result is illustrated in Fig. 3(a), and the vertical axis is for the output of the STN–GPe model represented by y3 (y3 = ySTN). In Fig. 3(a), solid black lines and dashed black lines represent stable and unstable equilibrium, respectively. HB is the Hopf bifurcation point, which can generate limit cycle oscillations (blue and red curves) corresponding to spontaneous beta oscillations in the model.

Fig. 3. (a) Codimension-one bifurcation versus C, where HB is Hopf bifurcation point. (b) Frequency versus C.

The bifurcation result in Fig. 3(a) shows that the system presents a unique stable fixed point (normal state) in the region of 0 < C < 6.056, and converges to a constant value for any initial condition, when C = 6.056, the model undergoes a Hopf bifurcation (HB1) point, and generates limit cycle oscillations. Thus it can be concluded that the STN–GPe model can be shifted from normal state to the abnormal spontaneous beta oscillations related to Parkinson disease when C is strong enough. It should be noted that in Fig. 3(a), the part of the curve where the model parameters are negative does not have any biological significance, but only serves as a mathematical description of the model.

Furthermore, the curve of frequency versus coupling strength C as illustrated in Fig. 3(b), suggests that increased coupling strength C can cause the oscillation frequency to increase. These results agree with those obtained by Liu et al.,[24] suggesting that simultaneously increased inter-coupling and intra-coupling strength can induce the beta oscillation to be generated and the oscillation frequency to be enhanced.

Next, we draw the diagrams for codimension-one bifurcation versus coupling strengths C1, C2, and C3, respectively, (and other parameters remain the default values as listed in Table 1) as shown in Figs. 4(a)–4(c). In Fig. 4(a), the bifurcation diagram of C1 shows that the system presents a bistable behavior (normal state) in the region of C1 < − 0.7706. In this parameter region, the output of the model converges toward one or the other stable fixed point, depending on the various initial conditions, and can switch between the two stable fixed points if perturbed. For −0.7706 < C1 < 5.812, the activity of the model corresponds to a unique stable fixed point. Finally, the model undergoes a Hopf bifurcation (HB) point when C1 = 5.812, and then generates limit cycle oscillations. In addition, the limit cycle oscillations in the red region may have a lower oscillation frequency than in blue region. This is because the larger the amplitude of limit cycle oscillation, the lower the oscillation frequency is. According to Figs. 4(a) and 4(b), we can conclude that the STN–GPe model can be shifted from normal state to the Parkinsonian oscillations state as long as the inter-coupling strengths C1 and C2 are strong enough. However, as shown in Fig. 4(c), increasing C3 can prevent the abnormal beta oscillations from being generated.

Fig. 4. Diagrams for codimension-one bifurcation versus (a) C1, (b) C2, and (c) C3, with other parameters remaining default values as listed in Table 1, BP representing branch point, and LP referring to limit point bifurcation. (d) Curves for frequency versus coupling strengths C1, C2, and C3.

The effects of coupling strengths C1, C2, and C3 on the frequency are shown in Fig. 4(d). These results suggest that the increasing of C1 and C2 reduces the frequency of oscillation, while the increasing of C3 facilitates the enhancement of the oscillation frequency. Thus we can conclude that compared with the inter-coupling strength between STN and GPe populations, the intra-coupling strength C3 of GPe population plays a contrary role in generating and regulating the spontaneous oscillations.

We further conduct theoretical analysis to explain why and how the increasing of C3 can enhance the generation of spontaneous neural oscillations and causes the oscillation frequency to increase.

According to Fig. 1(b), by approximating the sigmoid function in the GPe population with the first order Taylor expansion in the critical point

we can obtain the transfer function of the linearized model of the GPe population with inhibition self-feedback as follows:

where H(jω) is a second-order low-pass filter, and its natural frequency of GPe population with inhibition self-feedback is derived as

where is the natural frequency of GPe population without inhibition self-feedback (i.e., C3 = 0). Obviously, with the increase of C3, the natural frequency of the GPe population will become large, which contributes to the increase of spontaneous neural oscillations frequency.

Next, letting ω = 0 in Eq. (8), we can derive the equivalent gain of H(jω) as follows:

where γ′ is the equivalent gain of GPe population without inhibition self-feedback (i.e., C3 = 0), such that it can be concluded that the equivalent gain of GPe population becomes small with C3 increasing, which can prevent the model from generating the spontaneous oscillations.

3.1.2. Influence of interaction of coupling strengths on pathological beta oscillations

In this subsubsection, we explore the influence of the interaction of coupling strengths on pathological beta oscillations. Codimension-two bifurcation results versus coupling strength C1 and C2 are illustrated in Figs. 5(a) and 5(c) with different values of C3. As shown in Figs. 5(a) and 5(c), the behaviors of the model can be divided into three regions denoted as 1, 2, and 3 with different dynamic properties. The regions denoted as 1, 2, and 3 exhibit a bistable dynamic behavior, a unique stable fixed point and limit cycle oscillations state, respectively. For all parameter regions of C1 > 0 and C2 > 0, the STN–GPe circuit model undergoes the normal state corresponding to the region 2 and the beta band oscillations corresponding to the region 3. Hence, we mainly discuss the transition of dynamic states between region 2 and region 3.

Fig. 5. Influences of interaction of coupling strengths on dynamic behavior of the STN–GPe model, showing C2 versus C1 for C3 = 10 (a) and C3 = 30 (b), in which blue lines represent Hopf bifurcation curves, which correspond to HB point in codimension-one bifurcation diagram, and cyan lines are branch point curves, corresponding to BP point in codimension-one bifurcation diagram, and ((b) and (d)) frequency distribution diagrams.

In Figs. 5(a) and 5(c), the bifurcation results demonstrate that the increasing of inter-coupling strengths C1 and C2 can cause the abnormal beta oscillations to emerge. The boundary between stable region and unstable (limit cycle) region is in the form of a hyperbola, revealing that the product of C1 and C2 needs to be sufficiently strong for generating the beta oscillations. Figures 5(a) and 5(c) demonstrate that the increasing of intra-coupling (self-feedback) strength C3 results in the fact that the limit cycle oscillation regions in the two-dimensional plane of C1 and C2 become small, which is in accordance with Fig. 4(c).

To probe how the coupling strengths cooperate with each other to regulate the frequency of spontaneous oscillation, the frequency distribution diagrams of the spontaneous oscillations in the two-dimensional parameter space of C1 and C2 are obtained as shown in Figs. 5(b) and 5(d). Frequency distribution diagrams suggest that the spontaneous oscillations frequency decreases with C1 and C2 increasing. Besides, combining Fig. 5(b) with Fig. 5(d), we can draw the conclusion that the increase of C3 can increase the frequency of spontaneous oscillation in the beta oscillation region, which is consistent with Fig. 4(d).

It should be noted that the above results are in agreement with the results of Holgado et al.,[23] but inconsistent with those of Liu et al.,[24] where the increase of C1 results in the fact that the spontaneous neural oscillations in the STN–GPe model disappear (for the details see Fig. 4 in Ref. [24]). To verify the correctness of the present results shown as in Fig. 5, we further employ the describing function method to explore how the coupling strength exerts an effect on the generation of spontaneous oscillations.

The main idea of describing function method is to split a nonlinear system into two parts: nonlinear part N(A) and linear part L(jω). Solving the equation 1 + N(A)L(jω) = 0, i.e., L(jω) = −1/N(A), one can determine the limit cycle oscillations. For convenience, plotting the line of −1/N(A) and Nyquist diagram of L(j ω), if the two sides are crossed, there exists limit cycle oscillations in the nonlinear system.[43] As an effective method to analyze the limit cycle oscillations of nonlinear system, the describing function method is extensively used to gain an insight into the mechanism of neural spontaneous oscillations.[24,26]

In line with the study of Liu et al., we transform the sigmoid function in STN population into N(A), and the others in the STN–GPe model into linear part

with the sigmoid function in GPe population approximated by equivalent linear strength λ. We plot Nyquist diagram of L(j ω) and line of −1/N(A) in Fig. 6. The results demonstrate that the increase of coupling strengths C1 and C2 can enhance the generation of spontaneous in the STN–GPe model.

Fig. 6. Describing function analysis results of spontaneous oscillations with different values of C1 and C2, showing (a) line (yellow) of −1/N(a) and (b) Nyquist diagrams of L(jω).

Moreover, according to the principle of the describing function method, we can obtain the following function defining the stability boundary in the two-parameter plane (C1,C2):

Obviously, it is in the shape of hyperbolic curve, not a line as indicated by the results of Liu et al.,[24] which further demonstrates the correctness of the results as shown in Fig. 5.

In summary, the increasing of inter-coupling strengths C1 and C2, and the decreasing of intra-coupling strength C3 can cause the spontaneous beta oscillations to be generated in the STN–GPe model. It can be concluded that the STN–GPe circuit model can be shifted into a Parkinson-related oscillation state by changing the coupling strengths.

Previous studies demonstrated that increase of the coupling strength induces the the beta oscillation to be generated, as well as the oscillation frequency to be enhanced.[24] Our study reveals that the role of inter-coupling strength between STN and GPe populations played in generating the abnormal oscillations is different from that of intra-coupling strength of GPe, thus their influences on the oscillations should be separately studied. The present results reveal that the increase of inter-coupling strength between STN and GPe populations induces the spontaneous beta oscillations to be generated, and causes the oscillation frequency to decrease, while the increase of intra-coupling strength of GPe can prevent the model from generating the oscillations, and dramatically reduce the oscillation frequency. The present conclusions have some differences from those in previous study,[24] the main reason lies in the fact that Liu et al. combined inter-coupling strength and intra-coupling strength into one parameter (represented by C), thus drawing a confusing conclusion. Furthermore, we conduct a theoretical analysis to probe the influence of the inter-coupling strength of GPe population on the oscillations, revealing that the increase of inter-coupling strength leads the equivalent gain to decrease and the natural frequency to increase. The theoretical analysis results provide an insight into the generating mechanism and regulating mechanism of the oscillations induced by the inter-coupling strength of GPe population.

The study of Liu et al. claimed that the upper beta frequency in STN–GPe model originates from a high frequency cortical input, and the intra-coupling strength of GPe plays a transition role in determining low and upper frequency of the beta oscillations,[27] however the triggered mechanism is unclear. Our study reveal that the upper beta frequency can be generated by the model itself as long as the intra-coupling strength is large enough. Furthermore, the theoretical and simulation analysis further show that the intra-coupling strength can dramatically regulate the oscillation frequency by changing the natural frequency of the STN–GPe model. These results can shed light on the switching mechanism of the spontaneous beta oscillations triggered by the intra-coupling strength.

3.2. Influence of STN and GPe populations on pathological beta oscillations

In addition to the coupling strengths in the STN–GPe model, nonstructural changes in the STN–GPe circuit model, such as alterations of synaptic and intrinsic properties, are also critical for affecting dynamics behavior of the model.[27,44] In this subsection, we further explore how the STN and GPe populations exert influence on the spontaneous beta oscillations.

3.2.1. Influence of synaptic properties on pathological beta oscillations

First, we explore the influences of synaptic properties on the generation of pathological beta oscillations in the model. There are two parameters describing the synaptic properties of the STN–GPe model, i.e., excitatory and inhibitory average synaptic gains Hs and Hg. To explore how the two synaptic parameters interact with each other to induce abnormal Parkinsonian oscillatory behaviors, we keep the other parameter values listed in Table 1.

The bifurcation results are shown in Figs. 7(a) and 7(c) with different values of C3. In Fig. 7(a) with C3 = 10, the regions 1 and 2 represent the regions of normal state and the limit cycle oscillations, respectively. A Hopf bifurcation curve is the boundary to distinguish the two states. In Fig. 7(c) with C3 = 30, the regions 1 and 2 are normal states’ region, the region 3 represents beta rhythmic oscillations’ region, and the region 4 is the region of normal state and Parkinson’s state coexisting. The bifurcation results suggest that both increasing Hs and reducing Hg can induce the beta oscillations to be generated.

Fig. 7. Influence of interaction of synaptic gains on dynamic behaviors of the model, showing Hg with respect to Hs for C3 = 10 (a) and C3 = 30 (c). [(b), (d)] frequency distribution diagrams, [(e) and (f)] curves frequency with respect to Hs (setting Hg = 1) and Hg (setting Hs = 0.5) for C3 = 10 (e) and C3 = 30 (f). Other parameters remain default values as listed in Table 1. (g) Effect factor K of the parameters Hs and Hg on the generation of limit cycle oscillation. The black line in the codimension-two bifurcation diagram is the curve of the limit point bifurcation of cycles, corresponding to the LPC point in codimension-one bifurcation diagram.

As shown in Figs. 7(a) and 7(c), with intra-coupling strength C3 increasing, the beta oscillation region becomes small, which is in accordance with the scenario in Fig. 4(c). The effects of factors Hs and Hg on the generation of limit cycle oscillation are shown in Fig. 7(g), indicating that with the value of C3 small, the STN population plays a dominant role in generating spontaneous oscillation; while with the value of C3 large, the GPe population exerts a more important influence on producing the spontaneous oscillations. These results show that C3 plays a crucial role in reshaping the interacting relationship between STN and GPe populations.

To detect how Hs and Hg exert the influences on the frequency of beta oscillations, the frequency distribution diagrams in the (Hs, Hg) plane are given in Figs. 7(b) and 7(d), and the curves of frequency versus Hs (in the case of Hg = 1) and Hg (in the case of Hs = 0.5) with different C3 values are derived as shown in Figs. 7(e) and 7(f). The results indicate that both increasing Hs and reducing Hg cause the oscillation frequency to decrease. Therefore, it can be concluded that synaptic gain Hs and Hg in the STN–GPe model play key roles in triggering the beta band oscillations and regulating the oscillation frequency.

3.2.2. Influence of intrinsic properties on pathological beta oscillations

Finally, we argue that intrinsic properties, such as sigmoid parameters (i.e., es, eg, rs, and rg), are also capable of affecting the dynamic behaviors of the model, and induce pathological rhythmic activity related to Parkinson’s disease.

To verify our hypothesis, the diagrams of bifurcation and frequency distribution in both (es, eg) and (rs, rg) planes are derived to explore the influence of sigmoid elements on the generation of the spontaneous oscillations. Codimension-two bifurcation results with respect to es and eg, and with respect to rs and rg are illustrated in Figs. 8(a) and 8(c), and Figs. 9(a) and 9(c) with different values of C3, respectively. The bifurcation results suggest that there are four different oscillation regions 1, 2, 3, and 4, where region 1 represents bistable dynamic behavior region, regions 2 and 3 refer to a unique stable fixed point region and limit cycle oscillation region, respectively, and the region 4 (quite small) is the coexistence area of normal state and Parkinson’s state. On the whole, bifurcation results in Figs. 8 and 9 illustrate that the increase of es and rs enhance the generation of abnormal beta oscillations, while the decreasing of eg and rg cause the beta oscillations to be generated.

Fig. 8. Influence of maximum firing rate of the sigmoid function on dynamic behaviors of model, showing eg with respect to es for C3 = 10 (a) and C3 = 20 (c), [(b) and (d)] frequency distribution diagrams, eg with respect to es for C3 = 10 (e) and C3 = 20 (f). Other parameters remain default values as listed in Table 1, and (g) effect factor K of the parameters es and eg on generating limit cycle oscillation.
Fig. 9. Influence of slopes of sigmoid function on dynamic behaviors of the model, showing rg with respect to rs for C3 = 10 (a) and C3 = 20 (c), [(b) and (d)] frequency distribution diagrams, curves of frequency with respect to rs and rg for C3 = 10 (e) 10 and C3 = 20 (f). Other parameters remain the default values as listed in Table 1, and (g) effect factor K of the parameters rs and rg on generating limit cycle oscillation.

To quantitatively evaluate the role of intra-coupling C3, we further investigate the effect factor of sigmoid parameters on the generation of spontaneous oscillation as shown in Figs. 8(g) and 9(g). In the case of C3 = 10, both Kes and Krs are larger than 1; however, when C3 = 20, Keg and Krg are both larger than 1. These results suggest that with the value of small C3, the STN population plays a dominant role in generating the abnormal spontaneous oscillations, while with large value of C3, the GPe population exerts a more important influence on generating these oscillations.

The frequency distribution diagrams in the plane of (es, eg) and (rs, rg) are shown in Figs. 8(b) and 8(d), and Figs. 9(b) and 9(d), respectively. The frequency distribution diagrams reveal that the increasing of es and rs reduces the oscillation frequency, while the increasing of eg and rg can increase it. In Figs. 8(e) and 8(f), the curves of frequency with respect to es and eg for different C3 values show that both es and eg have important effects on the oscillation frequency, while in Figs. 9(e) and 9(f), the curves of frequency with respect to rs and rg with different C3 values demonstrate that rg plays a more crucial role in affecting the oscillation frequency than rs.

Previous studies showed that the STN primarily influences the generation of the beta oscillation while the GPe mainly determines its frequency. However, our study reveals that the intra-coupling strength of the GPe population plays a switching role in inducing the spontaneous beta oscillations through determining which population exerts a more important influence on the generation of the oscillations. For small value of the intra-coupling strength, the STN plays a dominant role in generating the beta oscillations; for large value of the intra-coupling strength, the GPe mainly determines the generation of the beta oscillations.

4. Discussion and conclusions

In the present study, we mainly employ the bifurcation analysis method to elucidate the mechanism of generating the pathological spontaneous beta oscillations triggered by coupling strengths and intrinsic properties of the STN and GPe populations. The results reveal that the increase of inter-coupling strength between STN and GPe populations enhances the generation of spontaneous beta oscillations, and results in the decrease of oscillation frequency. However, the increasing of intra-coupling strength of GPe can prevent the oscillations from being generated, and thus dramatically increasing the oscillation frequency. We further conduct theoretical analysis to provide an insight into the role of the inter-coupling strength of GPe population played in the generating and regulating of the oscillations. Furthermore, our study reveals that the intra-coupling strength of the GPe population can reshape the interacting relationship between STN and GPe populations, and determine which population exerts a more important influence on the generation of the spontaneous beta oscillations. The present results provide testable hypotheses for future experimental work.

In the present study, we focus on exploring the generation of the spontaneous beta oscillations caused by the excessive changes of parameters in the STN–GPe model. However, the inputs of the STN–GPe model from the cortex and striatum should also exert influence on the beta oscillations. In future study, we will further probe the roles of the inputs from cortex and striatum played in inducing physiological beta oscillations.

Reference
[1] Jankovic J 2008 J. Neurol Neurosurg Psychiatry 79 368
[2] Magnin M Morel A Jeanmonod D 2000 Neuroscience 96 549
[3] Thomas W Jesus S 2006 J. Neurophysiol. 95 2120
[4] Delong M Wichmann T 2007 Arch. Neurol. 64 20
[5] Marreiros A C Cagnan H Moran R J Friston K J Brown P 2013 Neuroimage 66 301
[6] Santens P Boon P Roost D V Caemaert J 2003 Acta Neurologica Belgica 103 129
[7] Hammond C Bergman H Brown P 2007 Trends Neurosci. 30 357
[8] Albada S J V Gray R T Drysdale P M Robinson P A 2009 J. Theor. Biol. 257 664
[9] Gillies A Willshaw D Li Z 2002 Proc. Biol. Sci. 269 545
[10] Holt A B Netoff T I 2014 J. Comput. Neurosci. 37 505
[11] Arvind K Stefano C Stefan R Ad A 2011 Front. System Neurosci. 5 86
[12] Mccarthy M M Moore-Kochlacs C Gu X Boyden E S Han X Kopell N 2011 Proc. Natl Acad. Sci. USA 108 11620
[13] Pasillas-Lepine William 2013 Biol. Cybern. 107 289
[14] Pavlides A John Hogan S Bogacz R 2012 Eur. J. Neurosci. 36 2229
[15] Avila I Parr-Brownlie L C Brazhnik E Castañeda E Bergstrom D A Walters J R 2010 Exp. Neurol. 221 307
[16] Moran R J Mallet N Litvak V Dolan R J Magill P J Friston K J Brown P 2011 Plos Comput. Biol. 7 e1002124
[17] Stein E Bar-Gad I 2013 Exp. Neurol. 245 52
[18] Tinkhauser G Pogosyan A Tan H Herz D M Kühn A A Brown P 2017 Brain Journal of Neurology 140 2968
[19] Chen M Guo D Li M Ma T Wu S Ma J Cui Y Xia Y Xu P Yao D 2015 Plos Comput. Biol. 11 e1004539
[20] Chen M Guo D Wang T Jing W Xia Y Xu P Luo C Valdes-Sosa P A Yao D 2014 Plos Comput. Biol. 10 e1003495
[21] Chen M Guo D Xia Y Yao D 2017 Front. Comput. Neurosci. 11 31
[22] Pavlides A Hogan S J Bogacz R 2015 Plos Comput. Biol. 11 e1004609
[23] Holgado A J N Terry J R Bogacz R 2010 J. Neurosci. 30 12340
[24] Liu F Wang J Liu C Li H Deng B Fietkiewicz C Loparo K A 2016 Chaos 26 123113
[25] Xia X F Wang J S 2014 Acta Phys. Sin. 63 140503 in Chinese
[26] Xu Y Zhang C H Niebur E Wang J S 2018 Chin. Phys. 27 048701
[27] Liu C Zhu Y Liu F Wang J Li H Deng B Fietkiewicz C Loparo K A 2017 Neural Networks: the Official Journal of the International Neural Network Society 88 65
[28] Grimbert F Faugeras O 2006 Neural Comput. 18 3052
[29] Touboul J Wendling F Chauvel P Faugeras O 2011 Neural Comput. 23 3232
[30] Spiegler A Kiebel S J Atay F M Knösche T R 2010 Neuroimage 52 1041
[31] Nevado-Holgado A J Marten F Richardson M P Terry J R 2012 Neuroimage 59 2374
[32] Jansen B H Rit V G 1995 Biol. Cybern. 73 357
[33] Mina F Benquet P Pasnicu A Biraben A Wendling F 2013 Front Comput. Neurosci. 7 1
[34] Wilson H R Cowan J D 1972 Biophys. J. 12 1
[35] Davidson C M de Paor A M Lowery M M 2014 IEEE Trans. Biomed. Eng. 61 957
[36] Gillies A Willshaw D 2007 Expert Rev. Med. Devices 4 663
[37] Hallworth N E Wilson C J Bevan M D 2003 J. Neurosci. Official J. Soc. Neurosci. 23 7525
[38] Hitoshi K Yoshihisa T Atsushi N Satomi C 2005 J. Neurosci. Official J. Soc. Neurosci. 25 8611
[39] Kita H 2007 Prog. Brain Res. 160 111
[40] Kita H Kitaia S T 1991 Brain Res. 564 296
[41] Kazanci F G Ermentrout B 2007 Siam Journal on Applied Mathematics 67 512
[42] Ashwin P Coombes S Nicks R 2016 J. Math. Neurosci. 6 2
[43] Phillips C L Parr J M 2012 Feedback Contron Systems 5 Beijing Science Press 482
[44] Mccarthy M M Ching S N Whittington M A Kopell N 2012 Curr. Opin. Neurobiology 22 693